/*
 * Let us call an integer sided triangle with sides a <= b <= c 
 * barely acute if the sides satisfy a^2 + b^2 = c^2 + 1.
 *
 * How many barely acute triangles are there with perimeter <= 25,000,000?
 */

#include <iostream>
#include "euler/gcd.hpp"
#include "euler.h"

BEGIN_PROBLEM(223, solve_problem_223)
	PROBLEM_TITLE("Almost right-angled triangles I")
	PROBLEM_ANSWER("61614848")
	PROBLEM_DIFFICULTY(2)
	PROBLEM_FUN_LEVEL(2)
	PROBLEM_TIME_COMPLEXITY("L*log(L)")
	PROBLEM_SPACE_COMPLEXITY("log(L)")
END_PROBLEM()

#if 0
// 2p(m+n)+1 <= L
static void test(int L, bool verbose)
{
	int trivial = (L-1)/2; // trivial solutions
	int count = 0;
	int dup = 0; // a,b,c>1, a,b,c all odd, a<>b.
	int hit = 0, miss = 0;
	for (int m = 2; m <= (L-1)/4; m++)
	{
		for (int p = 2; p <= (L-1)/(2*(m+1)); p++)
		{
			if (m%2==0 && p%2==0)
				continue;
			int q = 0, n = 0;
			if (euler::egcd(m, p, q, n) == 1)
			{
				++hit;
				n = -n;
				if (n < 0)
					n += m;
				if (q < 0)
					q += p;
				int a = m*p-n*q, b = m*q+n*p, c = m*p+n*q, l = 2*p*(m+n)+1;
				if (l <= L)
				{
					++count;
					if (m*q-n*p != 1)
						std::cout << "WRONG: ";
					if (verbose)
					{
						std::cout << "(m,n,p,q) = (" << m << "," << n 
							<< "," << p << "," << q << "), "
							<< "(a,b,c) = (" << a << "," << b << "," 
							<< c << "), a^2+b^2 = " << (a*a+b*b) << ", "
							<< "L=" << l << std::endl;
					}
					if (a!=b && a%2==1 && b%2==1)
						++dup;
				}
			}
			else
			{
				++miss;
			}
		}
	}
	std::cout << "# Enumerated pairs : " << (hit+miss) << std::endl;
	std::cout << "# Co-prime pairs   : " << hit << std::endl;
	std::cout << "# Perimeter < L    : " << count << std::endl;
	std::cout << "# Duplicate        : " << dup/2 << std::endl;
	std::cout << "# Non-trivial      : " << (count-dup/2) << std::endl;
	std::cout << "# Trivial          : " << trivial << std::endl;
	std::cout << "Answer             : " << (count-dup/2+trivial) << std::endl;
}
#endif

static int count_barely_acute_triangles(int L, int verbose)
{
	int trivial = (L-1)/2; // trivial solutions: a=1
	int count = 0;
	int dup = 0;

	// Callback function for solutions to mq - np = 1.
	auto f = [&](int m, int p, int q, int n) -> bool 
	{
		int l = 2*p*(m+n)+1;
		if (l > L)
			return false;
		if (m == 1 || p == 1) // skip trivial solutions
			return true;

		int a = m*p-n*q, b = m*q+n*p, c = m*p+n*q;
		if (verbose >= 2)
		{
			std::cout << "(m,n,p,q) = (" << m << "," << n 
				<< "," << p << "," << q << "), "
				<< "(a,b,c) = (" << a << "," << b << "," 
				<< c << "), a^2+b^2 = " << (a*a+b*b) << ", "
				<< "L=" << l << std::endl;
		}

		++count;

		// If a, b, c are all odd and a != b, then a^2+b^2 = b^2+a^2
		// correspond to different parameterizations of (m,n,p,q),
		// and one of them must be excluded.
		if (a!=b && a%2==1 && b%2==1)
			++dup;

		return true;
	};
	euler::generate_bezout_quadruples<int>(f);

	if (verbose)
	{
		std::cout << "# Trivial     : " << trivial << std::endl;
		std::cout << "# Non-trivial : " << (count-dup/2) << std::endl;
		std::cout << "Answer        : " << (count-dup/2+trivial) << std::endl;
	}
	return (count-dup/2+trivial);
}

// Todo: fix C4789 in divisor iterator.
static void solve_problem_223()
{
#if 0
	//const int L = 25; // 14
	//const int L = 50; // 32
	//const int L = 100; // 72
	const int L = 2500000; // 5352755
#else
	const int L = 25000000;
#endif

	int verbose = 0;
	std::cout << count_barely_acute_triangles(L, verbose) << std::endl;
}
